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Empirical mode decomposition is an adaptive signal processing method that when ap- 
plied to a broadband signal, such as that generated by turbulence, acts as a set of band-pass 
filters. This process was applied to data from time-resolved, particle image velocimetry 
measurements of subsonic jets prior to computing the second-order, two-point, space-time 
correlations from which turbulent phase velocities and length and time scales could be 
determined. The application of this method to large sets of simultaneous time histories is 
new. In this initial study, the results are relevant to acoustic analogy source models for 
jet noise prediction. The high frequency portion of the results could provide the turbulent 
values for subgrid scale models for noise that is missed in large-eddy simulations. The re- 
sults are also used to infer that the cross-correlations between different components of the 
decomposed signals at two points in space, neglected in this initial study, are important. 

I. Introduction 

Computational approaches to predicting jet noise may be categorized as acoustic analogy and numerical 
methods. The acoustic analogy methods use a rearrangement of the governing equations. The resulting 
equation contains a linear operator on one side of the equation that reduces to the wave equation for 
acoustic propagation at large distances from the source region. On the other side of the equation are terms 
that are significant within a relatively small region and are identified as the equivalent or analogous acoustic 
sources. The formal solution for the acoustic spectrum in the far field typically involves a convolution integral 
containing a wave propagator function and a correlation function of turbulence terms in the flow field . 1 
Numerical methods based on direct numerical simulation attempt to compute all the scales of turbulence in 
a flow followed by the computation of the associated radiated noise field. This type of computation requires 
a large amount of computer resources and time. Using less computer resources, large-eddy simulations 
compute the relatively larger scales of turbulence in the jet flow field. Consequently, higher frequency noise 
is missing from the resulting acoustic field spectrum calculations since the computational grid is too large 
to capture the turbulent noise sources at smaller subgrid scales. Thus, subgrid scale noise models have been 
proposed following the acoustic analogy approach . 2 Whether the full acoustic spectrum or just the high 
frequency portion of the acoustic spectrum is to be predicted, these methods require turbulence statistics 
from flow field measurements. Bodony & Lele , 3 using spatial filtering of a highly resolved, direct numerical 
simulation of a two-dimensional, low-Reynolds-number shear layer, computed the statistics for both the 
resolved, large-eddy-simulation-type scales and the unresolved, subgrid scales of turbulence. The parameters 
for the spatial filter had to be determined prior to its use. In this paper, empirical mode decomposition is 
used on measured data from jets to filter the data and separate the turbulent scales prior to the computation 
of turbulent statistics. The decomposition effectively filters the data in a manner that is automatic and 
signal dependent. 

Empirical mode decomposition (EMD) is a recently developed, adaptive signal processing method for any 
general, non-stationary, and nonlinear signal . 4 The method separates the data signal into a series of basis 
functions, called intrinsic mode functions, using the data itself. The intrinsic mode functions (IMFs) derived 
from signal decomposition have been used to distinguish physical phenomena of different frequencies and 
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wavelengths regardless of the long or short duration of the signal. For example, the method has been applied 
to study classical nonlinear systems, wind and water wave interactions, ocean waves and tides, tsunami 
waves, seismic waves, and atmospheric turbulence. 5 The newness of the method and its current issues, such 
as its limited mathematical foundation and uniqueness, 6 have not slowed its successful application. After 
a description of the velocity field measurements using time-resolved, particle image velocimetry (TR-PIV) 
and the jet operating conditions, a description is given of the EMD method as applied in this study. The 
application of EMD on a large set of simultaneous time histories generated by TR-PIV and subsequently 
computing two-point, space-time correlations is new. 

This paper provides some initial results determined from the second-order, two-point, space-time correla- 
tions and spectra computed using the total signals consisting of the flow velocity time histories from TR-PIV 
and the IMFs of the decomposed signals. Surveys of the axial fluctuating velocity signal and the IMF com- 
ponents of the decomposed signal space-time correlations are given at three reference points near the lipline 
of two subsonic jets. From the correlations, turbulence phase velocities and integral length and time scales 
are determined. For the higher speed jet of the two measured jets, the frequency dependent version of these 
turbulence values are presented and compared to the correlation derived values. It is the highest frequency 
IMF results for the phase velocities and the length and time scales that may be applicable to subgrid scale 
noise source modeling. Finally, as a prelude to future work, the importance of the cross-correlation between 
IMFs is considered. 


II. Description of Test Measurements 

An extensive set of measured jet noise and flow data has been acquired using the Small Hot Jet Acoustic 
Rig at the NASA Glenn Research Center. 7 Included in that set is data from time-resolved, particle image 
velocimetry (TR-PIV). Using this technique, velocity fields were measured in a jet flow at resolutions in 
both space and time suitable for computing spatial-temporal correlations and spectra. The details of the 
technique and the methods of data acquisition are found in Wernet 8 and in Bridges & Wernet. 9 Data from 
jets issuing from a converging nozzle with exit diameter D = 50.8 mm operating at exit Mach numbers of 
0.51 and 0.98 were used in this paper. The operating conditions for both jets, respectively labelled SP3 and 
SP7, are shown in Table 1 where T t /T, A, is the total temperature ratio relative to ambient conditions, Tg/Too 
is the static temperature ratio, Pt/Poo is the nozzle pressure ratio, Mj is the jet Mach number, Uj/coq is 
the jet acoustic Mach number, Uj is the jet exit velocity, and c a 0 is the ambient speed of sound. 


Case 

Tt/Too 

Ts/Too 

Pt/Poo 

M./ 

U J / Coo 

Uj{ m/s) 

SP3 

1.00 

0.95 

1.20 

0.51 

0.50 

172.8 

SP7 

1.00 

0.84 

1.85 

0.98 

0.90 

310.0 


Table 1. Test conditions for convergent nozzle flow measurements. 


As discussed in Wernet, 8 the size of the TR-PIV measurement field depends on the sampling or camera 
framing rate. The higher the sampling rate, the smaller the field of view becomes. With additional consid- 
erations such as the amount of storage space available for the data, the amount of data that can be obtained 
determines the length of the time histories for the TR-PIV measurements. The jet cases listed in the table 
were sampled at 25 kHz with a time history length of about 1 second resulting in the maximum Strouhal 
number for analysis of 3.67 for case SP3 and 2.05 for case SP7 with resolutions in Strouhal number of 0.014 
and 0.008, respectively. The measured field was 178.85 nun wide in the axial direction and 5.18 mm wide 
in the radial direction with a discretization of 70 by 5. The spatial resolutions are Ax/D = 0.0510 and 
Ar/D = 0.0255. Figure 1 shows a representation of the nozzle and the three TR-PIV measurement locations 
with equal scale in the r/D and x/D direction. 

Mean and turbulence quantities computed from the measured TR-PIV data are shown in Figures 2 
to 4. The mean axial velocity contours in Figure 2 reflect the wide disparity between the axial and radial 
dimensions. The negative sign on the radius is maintained by convention with the measurement coordinates, 
indicating that the measurements were made near the lipline below the centerline of a horizontal jet as shown 
in Figure 1. Line plots along each of the five radial measurement locations are shown for the normalized 
mean axial velocity in Figure 3. The lowest mean axial velocities occur along the line r/D = —0.49, nearest 
to the lipline of the jets. Figure 4 shows a similar set of line plots for the normalized axial turbulence 
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Figure 1. TR-PIV measurement domain size and location relative to the nozzle. 



(a) SP3 jet with Mj = 0.51, Ti/Too = 1.00 (b) SP7 jet with Mj = 0.98, Ti/7bo = 1.00 

Figure 2. Mean axial velocity contour plots. 




x/D x/D 


(a) SP3 jet with M j = 0.51, T t /T a 0 = 1.00 (b) SP7 jet with Mj = 0.98, T t /T x = 1.00 

Figure 3. Mean axial velocities at 5 equidistant radial locations from about r/D = — 0.39 to — 0 . 49 . Arrow points toward 
velocities at increasing radius. Dot marks reference point conditions. 




(a) SP3 jet with Mj = 0.51, Tt/T 0 0 = 1.00 (b) SP7 jet with Mj = 0.98, Tt/T a 0 = 1.00 

Figure 4. Axial turbulence intensities at 5 equidistant radial locations from about r/D — — 0.39 to — 0 . 49 . Arrow points 
toward turbulence intensities at increasing radius. Dot marks reference point conditions. 


intensities. Nearer to the nozzle exit the higher intensities occur near the lipline. These results are based 
on non-overlapping TR-PIV measurement locations. The middle location between about x/D of 6 and 
10 has mean and turbulence velocity measurements that appear to be consistent with those measured at 
the upstream location. The levels and slopes of the data lines at each radius appear connectable between 
measurement locations. The downstream measurement location is not consistent in the sense that the data 
lines are not all connectable across the gap between measurement locations. This was not resolved; hence the 
results from this location downstream should be considered approximate in terms of coordinate reference. 

Correlations and spectra for these data sets were computed for reference points nearest the lipline and 
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at the axial location of maximum difference in turbulence intensities across the shear layer and at the axial 
location nearest the end of the potential core. For the SP3 jet, these locations are at x/D = 2.39 and 
x/D = 6.90 and for the SP7 jet, x/D = 3.41 and x/D = 7.82. An additional reference point for each jet was 
at x/D = 10.0 in the downstream TR-PIV measurement location. The axial locations for these reference 
points are marked in Figure 3 for the mean axial velocities near the lipline and in Figure 4 for the axial 
turbulence intensities. 


III. Empirical Mode Decomposition 

Empirical mode decomposition (EMD) is an adaptive method to decompose any general non-stationary 
and nonlinear signal into a set of basis functions known as intrinsic mode functions (IMFs). 4 Flandrin et 
al. 10 describe the method as one of extracting the local oscillations from the local trend to obtain the IMF. 
The procedure is defined by an algorithm that is summarized as follows: 

1. Identify all maxima and minima of the signal. 

2. Interpolate between all maxima and all minima to define an envelope. 

3. Compute the average of the envelope at every point in the signal. 

4. Subtract the average from the signal to obtain the local oscillations. 

5. The local oscillation result is examined to see if it meets the following criteria for an IMF: 

(a) The difference between the number of extrema and the number of zero crossings is at most one. 

(b) The average of the envelope along the signal is zero or near zero according to some convergence 
criterion. 

If these conditions are not met, the current signal obtained at step 4 is used as the input signal starting 
at step 1 and the algorithm repeats until the criteria of step 5 are met. 

This process extracts one IMF C\{t) from the signal u(t ) leaving a residual Z\(t). The next IMF C^t) is 
then extracted from the residual Z\(t) and so forth until the residual contains one or no extreme value. The 
complete decomposition of the signal u(t) is achieved with a finite number of IMFs of the order N < log 2 K 
where K is the number of data points in the signal. The result is written as 

N 

u(t) = ^2c n (t) + z N (t) . (1) 

n = 1 

An example result of EMD performed on a TR-PIV time history is shown in Figure 5a. The total axial 
velocity (mean plus fluctuating) is shown as the top signal time history. Prior to performing EMD in this 
study, the total signal mean value was subtracted from the signal. The first process of removing the local 
oscillations from the local trend results in the first IMF containing the highest frequency content of the 
signal. As can be observed, each succeeding IMF contains lower frequency content as the local wavelength 
between zero crossings in the IMF increases. The residual is near zero since the total signal mean value 
was removed. Figure 5b shows the IMFs in greater detail over a fraction of the original signal length. The 
IMFs for a single time history are found to be nearly orthogonal 4 and uncorrelated. Further details and 
illustrations are given in Appendix A. 

The EMD method was developed in general for non-stationary and nonlinear signals and is often applied 
to signals with significant underlying trends or with short term events. Turbulence measurements such 
as those presented herein are broadband and approximately stationary in nature though they represent 
nonlinear behavior in the flow. Flandrin et al. 10 and Wu & Huang 11 discuss the application of EMD to 
broadband signals. The decomposition of a broadband signal into N IMFs gives a result that represents 
a signal having been processed by a dyadic filter bank with N filters. A dyadic filter bank is a set of 
overlapping, band-pass-type filters having a constant band-pass shape with each filter having half or double 
the frequency range of its neighboring filters. The mean frequencies of the filter bank are given by 

fn = fo2~ n ( 2 ) 
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(a) Complete measured time history and IMFs over 1 second. 
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(b) Expanded portion of time history over 0.06 seconds. 

Figure 5. Example TR.-PIV time history sampled at 25 KHz containing 24993 data points and all the corresponding 
intrinsic mode functions (IMFs). Data from the shear layer of the SP7 jet with Mj = 0.98, T t jT ^ — 1.00, at x/D = 3.41 
and r/D = —0.49. 


where f 0 is a constant and the number 2 is an approximate value. Figure 6a shows the power spectral 
density computed from the time histories shown in Figure 5 for the total fluctuating signal and each of the 
first seven IMFs. The spectra for the IMFs show the dyadic filter bank representation for the broadband 
data processed using EMD. Flandrin et al. 10 noted that there is no predetermination of the bandwidth and 
location of the filters, but they are automatic and signal dependent. However, for signals that have the same 
time history length and have similar broadband nature as those obtained from TR-PIV measurements, the 
filter bank will be similarly located for each signal in terms of the mean frequency and bandwidth of the 
filters. This is shown in Figures 6b and 6c at locations near the lipline further downstream. The original 
signal spectrum is recovered when adding the IMFs. 12 This is expected since from equation (1) the IMFs 
are added in the time domain prior to transforming to the frequency domain. As an illustration, Figure 6 
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shows the spectra computed after summing IMFs 2 and 3. The spectral level of the total fluctuating velocity 
is nearly recovered in this region of the spectrum. 




(a) x/D = 3.41 and r/D = —0.49 (b) x/D = 7.82 and r/D = —0.49 



(c) x/D = 10.0 and r/D = —0.49 

Figure 6. Power spectral densities computed for the axial velocity fluctuations at 3 axial locations near the jet lipline. 
SPT jet with Mj = 0.98 and Tt/Too = 1.00. Highest levels are the total fluctuating velocity spectra. Intrinsic mode 
function (IMF) spectra are labelled 1 to 7. Dashed lines are spectra for the sum of IMFs 2 and 3. 


IV. Correlation Definitions 


The data from TR-PIV is arranged as an array of time histories located at discrete points in a plane 
that cuts across the flow field of the jet. These data may be used to compute both second- and fourth-order, 
two-point correlations of the velocity fluctuations both for the total fluctuation and for the IMF components. 
However, the 1 second time histories were too short to obtain sufficient averages for reliable fourth-order 
correlations. Thus, this paper concentrates on computing second-order correlations. The second-order, 
two-point correlation is defined over a time length T by 


Rij(x,r 7 ,t) 


1 

r 


rT 


— / u' i {x,t)u'Ax + ‘n,t + T)At 


( 3 ) 


where u' (x, t) is the total fluctuating velocity (given that the fluctuating velocity is a sum of IMF components) 
in the i-th direction obtained from the total velocity using u[(x,t) = Ui(x,t ) — Ui(x,t), r/ is the spatial 
separation, r is the time delay, and an overbar denotes the time average of the quantity. (From hereon, 
the term ‘total’ refers to the total fluctuating velocity.) The normalized two-point correlation or correlation 
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coefficient is given by 


i(x,V,r) = ~ 


Rij(x,V,T) 

T Jo u i ( x > 0 T Jo" u j( x + V,t + T)u'j(x + r],t + T)dt 


1/2 


( 4 ) 


where the denominator integrals are the mean square values of the velocity fluctuations at the two locations 
in the field. 

Generalizing equation (1) for different fluctuating velocity components at each point in space 


N 


h (*£> G ^ (a?, f) T 2jAT (a?, f) , 


( 5 ) 


n=l 


we substitute this into equation (3) and multiply the terms to get the correlation equation that applies after 
empirical mode decomposition. 


Rij(x,ri,T) = — / Ci w (a;, t)C jn (x + rj, t + r) dt + ^ — / C iB (x, t) Cjfc(a; + r?, t + r) dt 

n=l ' n=l ' k=£n 

IV i r N l f T 

+ / z iN {x,t)Cj n {x + ‘n,t + r)dt + ^ — / Ci n (x, t)zjN(x + ri,t + t) dt 

n= 1 ' J ° n—1 ' J 0 

1 ^ 

+ — / z iN {x,t)z jN (x + rj,t + t) dt 
' Jo 


( 6 ) 


This equation is also normalized as in equation (4), thus all the normalized correlations involving the IMFs 
are relative to the root mean square values of the total fluctuating velocities at the two locations. 

To compute the length and time scales from the normalized, two-point correlation, we follow the approach 
of Kerherve et al. 13 The integral length scale is computed from 


and the integral time scale by 


r+oo 

A?j(®) = / rij(x,rik,T = 0) dr] k 

Jo 

r+oo 

Tij{x)= / nj(x,T1k = U ck T,T)dT 

Jo 


( 7 ) 


( 8 ) 


where U c k is the phase speed in the rjk direction. These scales are determined for signals that either contained 
all the frequency content of the total fluctuating velocity or the band-pass frequency content of the individual 
IMFs. To determine the frequency dependence of these scales, Kerherve et al. 13 used the complex coherence 
function. This function is computed using the Fourier transform with respect to the time delay of equation 
(3), the cross-power spectral density function, and the Fourier transforms of the velocity correlations at the 
two points in the field. See Kerherve et al. 13 for the details. The final results for the frequency dependent 
length scale and the frequency dependent time scale are 


r+oo 

A ij( x ’U)= dr] k 

Jo 


and 


r+oo 


Tij(x,u) = -- 

Uck{+>) 


\"fij(x,r]k,uj)\ drjk 


where, following Morris & Zaman, 14 


u c k(yj) =u/\d(j>(r) k ,u)/dr)k\, 


(9) 


(10) 


( 11 ) 


7 ij is the complex coherence function, 5ft denotes the real part, and (f>(r]k,i + ) is the phase of the complex 
coherence function. 


7 of 20 


American Institute of Aeronautics and Astronautics 



The approximately 1 second of data obtained from the TR-PIV measurements provided time histories 
containing 24993 samples. The calculation of the correlations and spectra were performed following pro- 
cedures given in Bendat & Piersol. 15 The time histories were divided into equal length segments. Each 
segment was windowed using a Kaiser-Bessel window (parameter a = 3.0), zero padded, and then processed 
by the fast Fourier transform. The segment transforms were summed and averaged to obtain auto- and 
cross-spectra. These were then inverse Fourier transformed to obtain the auto- and cross-correlations. With 
a fixed time history length, a trade-off has to be made between good frequency resolution and low variance. 16 
For this study, the time histories were divided into 194 segments with 256 points each with 50% overlap. 
This resulted in a frequency resolution of 48.8 Hz for spectrum estimates with a standard deviation of 8%. 

V. Results 

The velocity time histories at each point in the frame of TR-PIV data are decomposed into IMFs that 
are used to compute the second-order, two-point correlation of equation (6) normalized as in equation (4). 
The last three terms in equation (6) are correlations involving the residual of the decomposition. Figure 5 
shows the residual to be relatively small and nearly constant compared to most IMFs. Consequently, the 
integral involving residuals is small and the integrals involving residuals and IMFs are nearly zero. The 
latter follows from removing the nearly constant residual from the integral resulting in a computation of the 
average intrinsic mode value which is zero. Of the remaining two integrals, we will concentrate in this paper 
on computing the first integral; the second-order, two-point correlation of IMFs of the same mode number. 

Contour plots of ru(x, rji, r), the correlation between fluctuating axial velocity components, are shown 
in Figures 7 to 12 for the total fluctuating axial velocity and the first 5 IMFs from the two jets listed in 
Table 1 at the three reference points. The correlations are shown as a function of the normalized axial 
separation rji/D and the normalized delay time tUj/D. All the correlation plots are on the same contour 
scale from -0.1 to 1.0, thus the IMF correlation levels are relative to the total fluctuating velocity correlation 
level. A computed elliptic contour is included in the figures. These will be defined in the next paragraph. 
The IMF correlations all follow the same pattern. The highest frequency intrinsic mode IMF 1 is confined 
to a small region of space and time. As the frequency decreases with increasing IMF mode number, the 
correlation broadens in both space and time. There is also a noticeable change in the slope of the contours 
with IMF indicating a change in phase velocity with frequency. Two issues affect the results obtained from 
these correlations. The first is a lack of spatial grid points to properly resolve the rapid changes in the 
correlations in the high frequency IMF 1, especially in the lower velocity SP3 jet. The second is manifest 
as an anomaly visible in the contour plots at rft/D = 0 and centered on tUj/D = 0. It is especially visible 
in the generally lower contour levels in the higher number IMFs. It occurs within plus or minus one spatial 
grid point of rft/D = 0. Some results were affected by this as discussed below. 

To extract phase velocity and scale values from the correlation results, we fit an ellipse to a contour 
of the data, 3 shown as the black line shape in Figures 7 to 12. The contour chosen was 1/e times the 
correlation peak rn(:r,0,0) for the total fluctuating axial velocity or the individual IMFs. Assuming the 
correlation is monotonically decreasing, such as an exponential or Gaussian shaped function, an estimate for 
the integral length scale, L v = Ajq, following equation (7), was found. The distance is from the reference 
point, the origin in the (rji/D, tUj/D) plane, to the ellipse along ^\/D at tUj/D = 0. Similarly, the integral 
time scale, t v = r ii> equation (8), was estimated using the distance from the reference point to the ellipse 
along the line r/i/D = (U c /U j){tU j / D). The phase velocity U c /Uj was also determined from the ellipse 
equation. The details of this ellipse method are given in Appendix B. An advantage of this method is that 
it allows scale estimates to be made where measured data is lacking. For example, Figure 11a shows that 
the reference point is too close to the downstream edge of the measurement frame to allow the correlation 
to decay sufficiently along the constant phase velocity line to obtain an estimate for the integral time scale. 
However, with the ellipse equation determined from fitting the available data, an estimate for the time scale 
can be computed. The results using this method to estimate the phase velocities and the integral length 
and time scales are shown in Table 2 for the SP3 jet and in Table 3 for the SP7 jet. For reference, the 
mean axial velocities at the reference points are included in the tables. The values missing in the tables 
are due to the ellipse method giving unreliable results. Especially at the shear layer reference point and for 
higher number IMFs, the presence of the anomaly affects the location of the ellipse fitting points, skewing 
the tilt of the ellipse and affecting the length of the ellipse axes. Thus, the extracted values based on the 
ellipse method were inaccurate. As a check of the method in general, we compare the total fluctuating 
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x/D 


Total 

IMF 1 

IMF 2 

IMF 3 

IMF 4 

■u/Uj 

2.39 

U c /Uj 

0.705 

- 

- 

- 

- 

0.756 


LJD 

0.108 

- 

- 

- 

- 



Tr/Uj/D 

0.408 

- 

- 

- 

- 


6.90 

U c /Uj 

0.564 

0.709 

0.580 

0.550 

0.478 

0.594 


LJD 

0.388 

0.061 

0.137 

0.287 

0.509 



tJJj/D 

3.148 

0.135 

0.567 

1.686 

2.633 


10.0 

UJUj 

0.529 

0.623 

0.541 

0.498 

0.480 

0.548 


LJD 

0.459 

0.064 

0.144 

0.299 

0.529 



TrjUj/D 

3.638 

0.142 

0.655 

1.480 

3.089 



Table 2. Integral properties from correlations for the total fluctuating velocity and IMF components at three axial 
locations near the lipline r/D = —0.49 using the ellipse method. SP3 jet. Phase velocity, t/ c . Integral length scale, L v . 
Integral time scale, t v . Mean axial velocity, u. 


velocity results in Table 3 for the SP7, cold, Mj = 0.9 jet to those values found by Kerherve et al. 13 for an 
isothermal, Mj = 0.9 jet on the lipline near the end of the potential core. Using Table 3 at x/D = 7.82, 
we get U c = 185 m/s, L v = 18. 6mm, and t v = 0.50nrs. The Kerherve et al. 13 values are U c = 137m/s, 
L v = 19 mm, and t v = 0.62 ms. These results show that the ellipse method can produce correlation-based 
phase velocity and scale values comparable with other methods. The assumption here is that the correlation 
in this part of the flow follows an exponential or Gaussian shape implying little or no significant oscillations 
in the tail of the correlation. 

Within the limitations of the frequency range of the measurements and the nature of empirical mode 
decomposition, the IMF 1 values in Tables 2 and 3 represent the measured small-scale turbulent phase 
velocity and integral scales. Though the comparison is qualitative, since the measurements here are near 
the lipline of an axisymmetric jet, a couple of trends are similar to those found by Bodony & Lele 3 from the 
computed results for a two-dimensional shear layer. One is that the total turbulence phase velocities are less 
than the mean velocities on the higher-speed side of the shear layer and that the small-scale phase velocities 
may be higher (and in this case are higher) than the total phase velocities. The other is that the small-scale 
integral scales are much smaller than the total fluctuating velocity integral scales. The computed shear layer 
small-scale integral scales are 60 to 80% smaller in the downstream portion of the shear layer. For the two 
jet cases here, the small-scale integral scales are 80 to 95% smaller. 

Given that the intrinsic mode functions isolate a range of frequencies that decrease in frequency as the 
IMF mode number increases, the turbulent phase velocity and scale values listed in Tables 2 and 3 show the 
gross changes in these values with frequency. The phase velocity decreases with decreasing frequency and 


x/D 


Total 

IMF 1 

IMF 2 

IMF 3 

IMF 4 

u/Uj 

3.41 

UJUj 

0.687 

0.769 

0.540 

0.396 

- 

0.747 


LJD 

0.162 

0.090 

0.172 

0.252 

- 



tJJj/D 

0.984 

0.245 

0.529 

0.807 

- 


7.82 

UJUj 

0.596 

0.670 

0.586 

0.540 

0.384 

0.626 


LJD 

0.366 

0.076 

0.206 

0.375 

0.545 



tJJj/D 

3.057 

0.403 

1.036 

2.246 

2.266 


10.0 

UJUj 

0.568 

0.644 

0.562 

0.535 

0.460 

0.592 


LJD 

0.433 

0.090 

0.215 

0.401 

0.684 



tJJj/D 

3.606 

0.280 

1.089 

2.314 

3.476 



Table 3. Integral properties from correlations for the total fluctuating velocity and IMF components at three axial 
locations near the lipline r/D = —0.49 using the ellipse method. SP7 jet. Phase velocity, t/ c . Integral length scale, L v . 
Integral time scale, t v . Mean axial velocity, u. 
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Figure 7. Contours of m in the SP3 jet for the total fluctuating velocity and intrinsic mode functions (IMF) 1 to 5. 
Reference point x/D = 2.39 and r/D = —0.49. Black line is the computed ellipse. 



(d) IMF 3 (e) IMF 4 (f) IMF 5 

Figure 8. Same caption as Figure 7 with x/D = 6.90. 
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x Uj/D tU/D tU/D tUj/D 







Figure 10. Contours of rn in the SP7 jet for the total fluctuating velocity and intrinsic mode functions (IMF) 1 to 5. 
Reference point x/D = 3.41 and r/D = —0.49. Black line is the computed ellipse. 
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(d) IMF 3 


(e) IMF 4 


(f) IMF 5 


Figure 12. Same caption as Figure 10 with x/D = 10.0. 
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the integral length and time scales have smaller values at the high frequency (IMF 1) and increase in scale at 
progressively lower frequencies (higher IMF mode number). We next compare these results to those obtained 
for the frequency dependent length scale and the frequency dependent time scale following equations (9) and 
(10), respectively. 

The ability to compute the integrals in equations (9) and (10) without modeling the measured data is 
very difficult as noted by Kerherve et al. 13 In equation (9), the real part of the complex coherence decays 
slowly at low frequency with separation distance from the reference point and then begins to oscillate at 
large distances. As the frequency increases, the oscillations move closer to the reference point and may 
increase in intensity. Hence, within the limited spatial measurement frame, the spatial integration does not 
converge using only the measured data. Equation (10) uses the magnitude of the complex coherence which 
tends to only decay away from the reference point. Thus, above some frequency, the measurement space 
is large enough so that the integral should converge. Since the use of equation (9) with measured data is 
not possible, we use the decay of the complex coherence magnitude to estimate the frequency dependent 
length and time scales. Following Harper-Bourne 17 and Morris & Zaman, 14 the length scale is given by the 
distance at which the complex coherence magnitude decays by 1/e. (We acknowledge that 

Kerherve et al. 13 consider this to be a poor estimate.) For the time scale, we use the point at which the 
value \^ij (x, rji, <jS)\ / u c {uS) decays by 1/e. The frequency dependent phase velocity is computed directly from 
the phase of the complex coherence using equation (11). The results are shown in Figures 13 to 15 for the 
SP7 jet at the three axial locations near the jet lipline. 

Figure 13 shows the frequency dependent phase velocities. The results computed for the IMFs are shown 
in the peak region of the band-pass filter given approximately by the frequency range from /)) 2 -1 / 2 to 
2 1/2 , where /' is given by equation (2), except for IMF 1 which has the same upper frequency limit as 
the total fluctuating velocity. The phase velocities for IMF 1 and IMF 2 follow closely the phase velocity for 
the total fluctuating velocity except at the edges of the filter range. More scatter is found in the results for 
lower frequency IMFs 3 and 4 due in part to both insufficient frequency resolution and averaging. The color 
horizontal lines are the phase velocity values from Table 3 that were derived from the correlation results. 
The correlation derived phase velocities more accurately coincide with the frequency dependent results at 
the higher frequencies and at the downstream locations. The final comparison shows that the data follows 
the curve fit equation from Morris & Zaman 14 

^ = 0.062 In +0.701 (12) 

U.J \Uj J 

that has slightly different coefficients compared to the Harper-Bourne 1 ' version of this equation. The data 
from each of these locations in the SP7 jet would give slightly different fit slopes and offset values, but the 
basic form of the equation would be the same. 

The estimated frequency dependent length scales are shown in Figure 14. The basic trend is for the 
length scale to be somewhat constant at the low frequency followed by a decrease in length scale as frequency 
increases. These results follow those of Morris & Zaman 14 in that the Strouhal number where the length scale 





Figure 13. Phase velocity SPT jet with Mj = 0.98 and Tt/Too = 1.00 at 3 axial locations near the lipline r/D = —0.49. 
Symbols are frequency dependent phase velocities. Horizontal lines are correlation based phase velocities, eq. (16). 
Dashed line is eq. (12). 
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Figure 14. Length scale SP7 jet with Mj = 0.98 and Tt/Too = 1.00 at 3 axial locations near the lipline r/D = —0.49. 
Symbols are frequency dependent length scales. Horizontal lines are correlation based integral length scales, eq. (14). 
See Fig. 13c for legend. 


changes from constant to a decreasing value moves to a lower value as the reference point moves downstream. 
At the end of the potential core location shown in Figure 14b, there is insufficient axial extent to enable 
the computation of the lower Strouhal number length scales. The IMF results follow those for the total 
fluctuating velocity, though the estimates for the length scale are lower for IMFs 2 to 4. A possible reason 
for this will be discussed in the next section. A comparison of these results with the correlation integral 
length scale (horizontal lines for the values from Table 3) may indicate that the estimates for the frequency 
dependent length scale are too high. The oscillations inherent in the real part of the complex coherence 
would provide cancellations in the calculation of equation (9) resulting in smaller length scales that would 
coincide more closely with the integral length scales. 

Figure 15 shows the frequency dependent time scale for the total fluctuating velocity and the IMFs where 
there is sufficient data in the axial direction to make the calculation. The comparisons with the integral time 
scales are closer in value than found for the length scale comparison in Figure 14. Again, the IMF results 
follow the total fluctuating velocity results with more scatter in IMFs 2 to 4 especially in the shear layer. 
The accuracy of computing the integral time scale using the ellipse method was estimated for the IMF 1 
and IMF 2 results shown in Figure 15c. The 1/e point was found by interpolation along the constant phase 
velocity line of r n. The results are given by the dashed lines in the figure. For IMF 1, the directly computed 
integral time scale is 33% higher than the ellipse method time scale. With better resolution at IMF 2, the 
difference decreases to 11%. 





Figure 15. Time scale SP7 jet with Mj = 0.98 and T t /T = 1.00 at 3 axial locations near the lipline r/D = —0.49. 
Symbols are frequency dependent time scales. Horizontal lines are correlation based integral time scales, eq. (15). See 
Fig. 13c for legend. 
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VI. Discussion 


The results shown in the previous section are based on accepting that the IMFs of different mode numbers 
are nearly uncorrelated at a point and assuming that the neighboring point IMFs of different mode numbers 
are also uncorrelated. In this section, we will infer that the latter assumption is generally incorrect and that 
the second term on the right side of equation (6) cannot be ignored. This term represents all the two-point, 
cross-correlations in space and time between IMFs having different mode numbers. Figure 16a shows the 
individual two-point correlation as a function of the time delay at zero spatial separation for each IMF 
up to IMF 8 and the corresponding correlation for the total fluctuating velocity. The correlations are all 
normalized the same, as given in equation (4), to show the level of correlation among each individual IMF 
relative to the total fluctuating velocity correlation. Figure 16b shows incremental results from summing 
the correlations shown in Figure 16a. The IMF correlations sum to recover the total fluctuating velocity 
correlation at zero time delay which follows from the IMFs being orthogonal. However, when temporal delay 
occurs the individual IMF correlations no longer add to completely obtain the total fluctuating velocity 
correlation. The small difference between the total and the summed correlations indicates the the IMFs in 
this case are nearly uncorrelated at a point in space but not completely uncorrelated. Figure 17 shows that a 
small spatial separation between the two measurement points results in a much larger difference between the 
summed individual IMF correlations and the total fluctuating velocity correlation. The summed correlation 
reaches a peak that is 60% of the total correlation peak. Clearly, cross-correlations of IMFs with different 




Figure 16. Example m correlations as a function of the temporal separation with both measurement locations at the 
reference point. SPT jet. Reference point x/D = 7.82 and r/D = —0.49. 



Figure 17. Example m correlations as a function of the temporal separation with the measurement locations separated 
by r]i/D = 0.2041. SP7 jet. Reference point x/D = 7.82 and r/D = —0.49. 
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mode numbers are required to make up the difference. Further examples of correlations as a function of the 
axial spatial separation are shown in Figure 18 for zero time delay and in Figure 19 for a small time delay 
between measurement points. The requirement for including IMF cross-correlation terms is clearly evident 
as the sum of the individual IMF correlations falls short of the total fluctuating velocity correlation. 

The consequence of neglecting the cross-correlations between IMFs of different mode numbers can also 
be found in the phase velocity and scale results in Figures 13 to 15. For IMF 1, the effects appear to be 
small in the range of Strouhal number 1 to 2, the upper half of the measured frequency range. The total 
fluctuating velocity and IMF 1 results are in basic agreement. The IMF cross-correlation terms are not 
as important for the IMF 1 results except near the lower end of the frequency range. The IMF 2 to 4 
results occur in increasingly narrower frequency bands resulting in deviations in the phase velocities and 
decreases in the length and time scales compared to the total fluctuating velocity results. In Figure 14, 
the frequency dependent length scales of IMFs 2 to 4 are smaller than the corresponding length scales for 
the total fluctuating velocity. The difference comes from the neglected cross-correlation terms that pass 
through the linear process of Fourier transforming equation (6) and taking the real part in equation (9). The 
differences in the frequency dependent time scales are more difficult to discern since taking the magnitude 
of the complex coherence function, equation (10), is not a linear process. 




Figure 18. Example correlations as a function of the axial spatial separation with no time delay between measurement 

locations. SP7 jet. Reference point x/D = 7.82 and r/D = —0.49. 




Figure 19. Example m correlations as a function of the axial spatial separation with a time delay of rUj/D = 0.9764 
between measurement locations. SP7 jet. Reference point x/D — 7.82 and r/D = —0.49. 
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VII. Concluding Remarks 


A new application of the empirical mode decomposition (EMD) method has been shown in this paper. 
Instead of analyzing a single time history signal or a small set of time histories, a large array of closely 
spaced time histories derived from TR-PIV measurements have been analyzed. The results from EMD were 
equivalent to passing the time histories through a bank of band-pass filters. This data was then used in 
two-point correlation calculations resulting in the determination of turbulence phase velocities and length 
and time scales that only apply to a particular band-pass range. The highest frequency-range results may 
be applicable to subgrid scale noise source modeling. 

Some issues were identified that would need to be resolved to better use the method. 

• A higher sampling rate would increase the Strouhal number range that can be analyzed and/or allow 
higher speed jets to be studied. 

• The length of the time histories needs to be increased. This would allow for greater frequency resolution 
and/or allow a larger number of averages to be performed to reduce spectral variance and allow better 
results for fourth-order correlation calculations. 

• There needs to be finer spatial resolution in the TR-PIV measurements. The higher frequency corre- 
lations can change rapidly in a short distance. 

Unfortunately, hardware limitations may inhibit significant improvements in any of these issues. 

Finally, there is the issue of the importance of the cross-correlation of IMFs of different mode numbers. 
Since the IMFs represent data in overlapping frequency bands, it is not yet known if these space-time cross- 
correlations are between different frequency data or between the same frequencies in the overlapping regions. 
The fact that the highest frequency IMF 1 results appeared not to deviate much from the total results except 
at the low end of the frequency range suggests the latter explanation for the results. The lower frequency 
IMFs have overlap on both ends of the band-pass filter with results showing that the cross-correlation 
information is missing. Resolving these matters is part of future work. 

Appendix A — Some Details of Empirical Mode Decomposition 

Further details about the empirical mode decomposition (EMD) method as applied in this study are 
given in this appendix. Specifically, information is given about interpolation and sifting, stopping criteria, 
and boundary conditions that were used in the current work. Any issues related to these are discussed in 
the references cited below 

Empirical Mode Decomposition 

The EMD method may be applied to any general, non-stationary, and nonlinear signal. 4 The signal u(t) 
is sampled to create a set u(i) where i = 1,2, ... ,K. We start by identifying the extrema of u(i), the 
positive peaks and the negative peaks in our case since the mean value of the signal is removed from the 
signal prior to decomposition. The positive peaks are used to construct an interpolating function e maa; (i) 
using a cubic spline. Similarly, the negative peaks are used to construct the interpolating function e rnin (i). 
The original signal u{i) and e max (i) and e m j n {i) are illustrated in the top part of Figure 20. As can be 
seen, the interpolating functions represent the envelope of the signal u(i). We now compute the average of 
the envelope m(i) = {e max (i) + e m i n (i))/2, which gives the local trend in the signal. Subtracting the local 
trend from the signal results in the local oscillation. Ideally, the local oscillation is a simple function with 
varying amplitude and frequency with zero mean value and with the number of extrema and the number of 
zero crossings differing by no more than one, called an intrinsic mode function (IMF). However in practice, 
these conditions do not immediately occur and a process of iteration called sifting is followed. The current 
local oscillation becomes the new signal for which new envelope functions are determined, labelled ‘Sift 1’ in 
Figure 20, and the average of the envelope is computed. This new local trend is subtracted from the current 
signal to obtain the new local oscillation. As noted in the figure, at least 8 siftings are required for the local 
trend to approach zero and the resulting IMF C(i) to have zero mean. It should also be noted that the 
envelope of the IMF is approximately symmetric about the axis. 
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Figure 20. Example of the process of empirical mode decomposition. 


Stopping Criteria 

As stated in the main text, the criteria that a local oscillation must meet for it to be considered an IMF are: 

1. The difference between the number of extrema and the number of zero crossings is at most one. 

2. The average of the envelope along the signal is zero or near zero according to some convergence criterion. 

To achieve these conditions in the numerical calculations, we followed the stopping criteria of Rilling et al. 18 
First of all, the sifting will stop any time the number of extrema is less than or equal to 2. With sufficient 
extrema in the current local oscillation, we compute a relative mean amplitude at each data point. 

/ .\ €-max{i) T 

= p— 7W 

If the denominator is equal to zero, set s(i) = 0. To stop sifting, both the following must be true: 

1. s(i) < (f> 2 at all points. 

2. Let k be the number of points where s(i) < i, then k/K > 1 — a. If a = 0.05, then 95% of the points 
meet the condition s(i ) < (f> 

To achieve these conditions, we must have <f >2 > 4>i- The recommended values of </> 2 = 0.5 and <f>\ = 0.05 
given by Rilling et al. 18 were used. 

Even though a stop-sifting criterion is specified, sifting can go on for many iterations potentially leading 
to undesirable results. Wu & Huang 19 discuss the issues involved with the inexactness of computing a 
stopping criterion and suggest setting the number of sifts to a fixed value of 10. We followed this by limiting 
the number of sifts to 10, but accepting less sifts if the above enumerated criteria were met. 

Boundary Conditions 

To obtain the envelopes of the signals shown on Figure 20, boundary values must be specified at the ends of 
the signal since typically these are not extreme points. We followed the procedure to extrapolate from the 
two extrema of the same type nearest the end point. 19 For the maxima envelope, if the extrapolated value 
is higher than the end point, choose the extrapolated value. Otherwise, keep the end point value. Similarly 
for the minima envelope, use the extrapolated value if it is lower than the end point value. 

18 of 20 


American Institute of Aeronautics and Astronautics 


Appendix B — Statistics from the Ellipse Equation 


The contours of a two-point, second-order correlation of the axial fluctuation velocity has an elliptic 
shape in the (771/ D,tU j / D) plane, see Figure 8a for an example. To determine an equation for an ellipse, 
we extract coordinates for the points in the plane where the level is 1/e times the peak of the correlation 
using interpolation on the data grid. These points are fit to the ellipse equation 


(V 1 N 

) 2 + b ^ tUj +cl 

( tUj\ 

\D , 

1 + D D +C ' 

v D J 


+^ + 


D 


+ f = 0 


(13) 


using the least-squares, non-iterative approach from Halff & Flusser 20 to obtain the coefficients a through 
/. Note from the diagram in Figure 21 that the center of the ellipse ( x c ,y c ) does not necessarily coincide 
with the origin of the axis system. 

Given the coefficients, we estimate the length scale as the distance from the origin to the ellipse along 
the line tUj/D = 0. Substituting this into equation (13) and solving for r/i/D results in 


L 

D 




(14) 


The time scale is found from the distance from the origin to the ellipse, where the semi-major axis a x 
intersects the ellipse, and projected onto the tUj/D axis. Figure 21 shows this distance to be 


tUj 

D 


= a x sin a + y c 


The phase speed is determined from the angle of rotation the ellipse makes about the origin 


1 


U r . 


Uj tan 9 


(15) 


(16) 


Using analytic geometry, equations are derived to aid in the solution of equations (15) and (16) based on 
the coefficients in equation (13). 

a x sin a + y c 


tan0 = 

a x cos a + x c 
— tan 2 a) 

a x = \ — 5 

y a — c tan a 
/i = ax 2 c + bx c y c + cy 2 c - f 


x c = 


Vc = 


2 cd—be 
b 2 — 4 a c 
2 a e — b d 
b' 2 — 4 a c 


tan a = — 


a — c 

b 


a — c 

b 


+ 1 


(17) 

(18) 

(19) 

(20) 

(21) 

(22) 
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Figure 21. Diagram of ellipse system for computing phase velocity and length and time scales. 
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